Isovector splitting of nucleon effective masses, ab-initio benchmarks and extended 
stability criteria for Skyrme energy functionals 
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We study the effect of the splitting of neutron and proton effective masses with isospin asymmetry 
on the properties of the Skyrme energy density functional. We discuss the ability of the latter to 
predict observable of infinite matter and finite nuclei, paying particular attention to controlling the 
agreement with ab-initio predictions of the spin-isospin content of the nuclear equation of state, as 
well as diagnosing the onset of finite size instabilities, which we find to be of critical importance. 
We show that these various constraints cannot be simultaneously fulfilled by the standard Skyrme 
force, calling at least for an extension of its P-wave part. 

PACS numbers: 21.30.Fe, 21.60.Jz 



I. INTRODUCTION 

In the study of medium to heavy mass nuclei, nuclear 
Energy Density Functional (EDF) approaches, based on 
self-consistent Hartree-Fock (HF) methods and their ex- 
tensions, constitute the theoretical tool of choice 0. 
Thanks to the development of better energy functionals 
and to the increase of computer resources, nuclear EDF is 
on the edge of becoming a predictive theory for all nuclei 
but the lightest. This is not only true for ground state 
properties, such as binding energies, radii or multipoles 
of the density, but also for low-energy spectroscopy and 
decay probabilities 0. 

However, the accuracy and predictive power needed 
for unknown regions of the nuclear chart still leave a lot 
of room for improvement. The phenomenological nature 
of Skyrme functionals makes their ability to faithfully 
predict observable or phenomena not linked with those 
used for their construction quite weak. Indeed, the lim- 
ited number of adjustable parameters (compared to the 
wealth of nuclear observable to be matched) turns fit- 
ting a Skyrme functional into an overconstrained prob- 
lem (which, of course, does not prevent some parts of it 
from being underconstrained) . 

As a direct consequence, many properties of existing 
parameterizations are biased to the fitting procedure and 
the limited analytical form of the Skyrme force, rather 
than to physical reasoning. A well-known example is the 
equation of state (EOS) of Pure Neutron Matter (PNM), 
which is sometimes subject to a pathological collapse at 
high density when not explicitly constrained. This is 
problematic insofar as one of the major challenges of con- 
temporary nuclear theory is to predict properties of very 
isospin asymmetric nuclear systems, i.e. neutron rich nu- 
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clei and matter in neutron stars. Experimental data be- 
ing unavailable in this domain of isospin, one has started 
relying on ab-initio theoretical results to constrain isovec- 
tor properties of the functional. It has led to the con- 
struction of the "Saclay-Lyon" SLy series of parameter- 
izations 0, H || by fitting (among other quantities) a 
theoretical equation of state of neutron matter. 

Isovector features of the nuclear EOS are crucial for 
a good understanding of neutron stars, exotic nuclear 
collisions produced at radioactive beam facilities and to 
describe the structure of exotic nuclei. For instance, the 
density dependence of the volume symmetry energy de- 
termines the proton fraction in f3 equilibrium in neutron 
stars, which ultimately drives the cooling rate and neu- 
trino emission |5| . The high-density part of the symmetry 
energy, which happens to be strongly model dependent, 
also influences significantly the isospin diffusion in heavy- 
ion collisions . Finally, the low-density part of the sym- 
metry energy is correlated with the size of neutron skins 
in finite nuclei . 

Beyond global isospin-dependent properties of the 
EOS, the isovector part of nucleon-dependent quanti- 
ties may influence the behavior of the above mentioned 
systems. Thus, collision observable depend on the mo- 
mentum dependence of the mean-field, in particular 
on its isovector component @, Also, some prop- 
erties of neutron stars require a precise know ledg e of 
isoscalar and isovector nucleon effective masses [T(J, [H| ■ 
The latter, which drives the splitting of neutron and 
proton effective masses with neutron/proton asymme- 
try, will serve as a starting point for our study. In- 
deed, a lot of efforts has recently been devoted to the 
microscopic characterization of neutron and proton ef- 
fective masses in infinite Asymmetric Nuclear Matter 
(ANM) IS [13, Q, [IE EE E3 EE EM 123, Either in 
ANM or in nuclei, the two species acquire different effec- 
tive masses. This property is quantified by the difference 
Am* (J) =m*(I)-m*(I), where I = (p n -p P )/(Pn + P P ) 
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is the isospin asymmetry while p n and p p denote neutron 
and proton densities, respectively. Note that the differ- 
ent effective masses m* discussed in the text always refer 
in fact to the ratio m*/m, where m is the bare nucleon 
mass. The latter is taken to be the same for neutrons 
and protons. 

This effective-mass splitting, though, is only one of a 
wealth of observable which can be subject to comparison 
between ab-initio predictions and EDF models. In this 
work we present results of a classical yet long unused 
test: the separation of infinite Symmetric Nuclear Mat- 
ter (SNM) potential energy per particle into spin-isospin 
channels. 

We shall also pay particular attention to controlling 
instabilities (i.e. non-physical spontaneous breaking of 
spin, isospin and/or spatial symmetries), and correlate 
Am* (I) with vector properties of the functional. We 
thus investigate the behavior of the latter with respect 
to the breaking of time-reversal invariance and the onset 
of spin polarization, looking for an overall consistency 
check of its spin-isospin content. Indeed, such proper- 
ties will become more and more important as one at- 
tempts to use full-fledged Skyrme functionals to study 
odd-mass nuclei, calculate rotational properties through 
self-consistent cranking calculations, or use more general 
dynamical methods |22j . 

This paper is organized as follows: in section [H] we 
present the set of Skyrme parametcrizations used in this 
study and examine basic properties of nuclear matter and 
finite nuclei. From then, in section II I II we perform a 
more detailed study of the spin-isospin content of the 
functionals and of their stability against finite-size spin 
and isospin perturbations using response functions in the 
Random Phase Approximation (RPA). 



II. CONSTRAINING THE ISOVECTOR 
EFFECTIVE MASS 

The nucleon effective mass to* is a key property char- 
acterizing the pr opag ation of (quasi)nucleons through the 
nuclear medium [23j . It is a reminder of the non-locality 
and energy dependence of the nucleon self-energy S(fe, e), 
themselves originating from the finite range and non- 
locality in time and space of the nucleon-nucleon inter- 
action. Mean-field-like theories of finite nuclei or infinite 
matter assume an on-shell propagation of the nucleons. 
The effective mass associated with such an on-shell prop- 
agation does not take into account the fragmentation of 
the single-particle strength and thus, a limited part of the 
effects associated with the energy dependence of e). 
Finally, the total energy is calculated by considering the 
quasi-holes (particles) to have spectroscopic factors of 1. 
In this context, either microscopic [24| or making use of 
phenomenological interactions or functionals mean- 
field methods do not correspond to a naive Hartree-Fock 
theory and always amount to renormalizing a certain 
class of correlations into the effective vertex. However, 



the energy dependence of the self-energy arising from 
the correlations only influences the position of the quasi- 
particle peak energy. The standard nuclear EDF meth- 
ods differ in several respects from Kohn-Sham Density 
Functional Theory [25|. In particular, the strategy usu- 
ally followed in the former is to leave room for improve- 
ment through further inclusions of correlations. This can 
be done by performing (Quasiparticle)-Random-Phase- 
Approximation (QRPA) calculations [27J or by em- 
ploying the Generator Coordinate Method (GCM) on top 
of symmetry-restored mean- field states j2g • 

Thus, the effective mass adjusted at the pure mean- 
field level is not expected to generate single-particle spec- 
tra matching exactly experimental data extracted from 
neighboring odd-mass nuclei. In particular, the coupling 
of single-particle motions to surface vibrations in closed- 
shell nuclei is known to increase the density of states at 
the Fermi surface and thus the effective mass [2^, [3(]] • A 
mean-field isoscalar effective mass to* lying in the interval 
0.7/0.8 in SNM, is able to account for a good reproduc- 
tion of both isoscalar quadrupole giant resonances data 
in doubly closed-shell nuclei [3l| and of single-particle 
spectra in neighboring ones provided particle-vibration 
coupling has been properly included. When the latter 
coupling is taken into account, the effective mass be- 
comes greater than one for states near the Fermi surface. 
Certainly, a lot remains to be done to understand these 
features microscopically in more involved cases |32^ . This 
is not only true for mid-shell nuclei where the coupling 
to both rotational and vibrational states can be impor- 
tant, but also for exotic nuclei where the coupling to the 
continuum becomes crucial and where shape coexistence 
and/or large amplitude motion appear more systemati- 
cally. 

In very exotic systems, the isovector behavior of to* 
and to* should play an important role. However, so far, 
no experimental data from finite nuclei has allowed a de- 
termination of the effective mass splitting as a function of 
neutron richness. In this context, ab-initio calculations 
of ANM are of great help. Non-relativistic Brueckner- 
Hartree-Fock (BHF) calculations, with or without three- 
body force, and, with or without rearrangement terms 
in the self-energy, predicted Am* (I) to be such that 
m n > m p in neutron-rich matter, that is, for J > 0. 
Such a conclusion was also reached by calculating the 
energy dependence of the symmetry potential (the Lane 
potential |33|) within a phenomenological formalism 
The latter result was confirmed by microscopic Dirac- 
Brueckner-Hartree-Fock (DBHF) calculations [34|. The 
situation regarding the prediction of the effective mass 
splitting was complexified due to an apparent contra- 
diction between results obtained from BHF 0, [w| and 
DBHF calculations Ha. However, the situation was fi- 
nally clarified in Ref. |l9l l20j where the importance of 
the energy dependence of the self-energy and the need to 
compare the non-relativistic effective mass with the vec- 
tor effective mass in the relativistic framework [35j | were 
pointed out. 
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Thus, the sign of the splitting is rather solidly pre- 
dicted. However, its amplitude is subject to a much 
greater uncertainty. Starting from that observation, the 
goal of the present section is to study the impact of the 
effective-mass splitting on properties of exotic nuclei pre- 
dicted by Skyrme-Hartree-Fock calculations. As far as 
the effective-mass splitting is concerned, one expects con- 
sequences onto structure properties of neutron-rich nu- 
clei. As a relatively large asymmetry may be necessary 
to reveal the influence of the splitting, data from nuclei 
not yet studied experimentally should provide crucial in- 
formation in that respect. As the effective mass gov- 
erns the density of states at the Fermi surface (together 
with the spin-orbit and the tensor forces), the amplitude 
of the splitting may influence properties such as masses 
and single particle properties of exotic nuclei, the evolu- 
tion of isotopic shifts across neutron rich closed-shell nu- 
clei or shell corrections in superheavy nuclei around the 
(N = 184, Z = 120) island of stability |M El S HI ■ 
Also, neutron and proton correlations beyond the mean- 
field should develop rather differently depending on the 
direction and amplitude of the effective-mass splitting. 
This could be true for static and dynamical pairing corre- 
lations as well as for the coupling to vibrational and rota- 
tional states. Finally, the effective mass splitting should 
leave its fingerprint onto the characteristics of isovector 
vibrational states of different sorts in neutron rich nu- 
clei EH. 



A. Fitting protocol 

Trying to keep a coherence in the way we construct 
Skyrme functionals, we take the fitting protocol used to 
define the SLy forces Q as a basis for the present 
work. Also, we pay attention to the fact that any im- 
proved or complexified functional includes all features 
validated by the SLy ones. 

The SLy functionals were derived from an effective 
interaction, that is, time-odd components of the func- 
tional are eventually linked to time-even ones. How- 
ever, some terms in the functional given by Eas. lA6alA6d1 
were dropped for some of the parameterizations. For in- 
stance, time-odd terms of the form s q ■ As g / have not 
been considered when calculating odd nuclei or ro- 
tational states |4^. Also, J 2 contributions to the energy 
functional associated with momentum-dependent terms 
in the central Skyrme force have been omitted for some 
of the SLy parameterizations, as it has been the case for 
most of the Skyrme functionals so far [l|. Rigorously, 
omitting terms is inconsistent with the idea of deriving a 
functional from an interaction. In any case, the latter ap- 
proach is only used for simplicity until proper adjustment 
or derivation of time-odd terms is feasible. We refer the 
reader to Refs. for a more extensive discussion 

on the subject. We presently take the SLy5 parameter- 
ization as a starting point. Thus, the two-body part of 
the center of mass correction is omitted whereas the J 2 



terms are fully kept. The spin-orbit term is the standard 
one, with a single parameter adjusted on the splitting of 
the 3p neutron level in 208 Pb. 

Within this general scheme, we have built a series of 
three new Skyrme forces, denoted hereafter /— ,/o and 
/+. The departures from the SLy protocol considered 
presently are (i) a better control of spin-isospin instabil- 
ities via Landau parameters (ii) the use of two density- 
dependent zero-range terms |43j (iii) a constraint on the 
isovector effective mass, such that, in neutron-rich sys- 



for /_ 



for fo and m* > 



terns, m* < 
for /+. 

With two density dependent terms, the compressibil- 
ity and the isoscalar effective mass are no longer bound 
together and can be chosen independently. However, 
this is not directly used here and an isoscalar effective 
mass of m* = 0.7, close to the SLy5 value, is chosen for 
the three parameterizations /_,/o,/+. The additional 
freedom brought about by the second density-dependent 
term is only used to adjust more easily the high-density 
part of the PNM EOS (see below). In the end, the only 
parameter subject to variation between /_, /o and /+ is 
the isovector effective mass m* which, m* being constant, 
drives the splitting Am* (/) . 

In the present work, we use the SLy5 force as a ref- 
erence, and include a comparison with the recently pro- 
posed LNS parameterization [44| which was also built to 
match the splitting of effective masses and the neutron 
matter EOS predicted by BHF calculations. The SkP 
force initially built for the study of pairing effects, 
will be used for a special purpose in the discussion about 
instabilities. 



B. Elementary properties of studied forces 



As we focus on the behavior of effective masses m* with 
isospin asymmetry, we recall that these quantities are re- 
lated to the dependence of the energy density functional, 
Eas. IA6allA6dl on kinetic densities r g , as 

h 2 

+ C T a po + ql CI po (la) 



2m* (I) 



m 



&H 

dr„ 



2m 



m*(I) 



m mm 

— + qi — * 



(lb) 



where po is the scalar-isoscalar density and q = +1,-1 
respectively for neutrons and protons (for the definition 
of CI coefficients, see appendix IA"|) . The splitting of ef- 
fective masses, quantified by 



Am* (I) _ m* (/) m* p (I) 



(2) 



m mm 
is governed by the isoscalar and isovector effective masses 
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771* 
771 

771* 



i + 1? (C5 - ci) p 



1 + K s 
1 + Kv 



(3a) 
(3b) 
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We use the usual convention for the isovector effective 
mass, which stems from its definition through the en- 
hancement factor k v of the Thomas-Reiche-Kuhn sum 
rule |46|. However, m* and n v are not isovector quanti- 
ties in the sense of isovector couplings of the functional. 

In the following, we shall discuss the value of Am* (7) 
at I = 1, which we note Am* in the following, for the 
sake of brevity. We have 

Am* _ 2(« v - k s ) 

m (1 + k s ) 2 — (k v — k s ) 2 ' 

such that Am* > for n v > k s , or equivalently m* < m*, 
or CI < 0. 



Table I: Infinite nuclear matter properties of the Skyrme 
forces quoted in the text. The quantities p sat and E/A de- 
note the density and energy per particle at saturation in SNM. 
The symmetry energy and the compressibility (for symmetric 
matter) are respectively 32 MeV and 230 MeV for SLy5 and 
all f x forces. In the case where m* ~ 0.7, k s ~ 0.43, so we 
have Am* > if k v > 0.43. 



Force 


Psat 


E/A 


m* 


ft V 


m*. 


Am* 


SLy5 


0.161 


-15.987 


0.697 


0.25 


0.800 


-0.182 


f- 


0.162 


-16.029 


0.700 


0.15 


0.870 


-0.284 


fa 


0.162 


-16.035 


0.700 


0.43 


0.700 


0.001 


f+ 


0.162 


-16.036 


0.700 


0.60 


0.625 


0.170 


LNS 


0.175 


-15.320 


0.825 


0.38 


0.727 


0.227 


SkP 


0.170 


-16.590 


1.030 


0.32 


0.760 


0.418 



Bulk properties of f x parameterizations are displayed 
in Table Q] We note that, while the position of the sat- 
uration point varies little between our forces (SLy5 and 
f x ), this consistency is lost in the case of LNS and SkP. 
These properties depend on the observable used in the fit- 
ting procedure. In the case of LNS, the saturation point 
relates to an Extended Brueckner-Hartree-Fock (EBHF) 
calculation 0], predicting values of (E/A) sat and p sa .t 
which are larger than empirical ones. A similar but lesser 
trend is observed for SkP. In this case it seems to be cor- 
related with the choice of effective masses and their inter- 
play with other parameters of the force. Indeed, binding 
energies computed with SkP compare satisfactorily with 
experimental ones, while LNS suffers in this respect from 
the lack of readjustment of the saturation poin t on nu- 
clear data. As it has been shown in Ref. [43, nuclear 
binding energies are highly sensitive to the choice of the 
energy at saturation, which is therefore constrained to a 
very tight interval if one wants to reproduce such quanti- 
ties. This constraint is especially tight compared to the 
uncertainty of ab-initio predictions. Despite the fit of 
surface properties (C^ p parameter) on a set of nuclear 
data, the accuracy of binding energies predicted by LNS 
is of the order of 5%, to be compared with less than 1% 
for SLy5. 



C. Properties of the nuclear matter EOS 

It is interesting to note that SLy parameterizations 
were fitted to PNM EOS with the idea of improving 
isospin properties of the functionals. One consequence 
was to generate functionals with Am* < 0, in opposition 
to ab-initio predictions. On the other hand, older func- 
tionals such as SIII [48| and SkM* [43, which were not 
fitted to PNM, had Am* > 0. The same exact situation 
happens for the Gogny force [50T |. Thus, improving global 
isovector properties (EOS) seems to deteriorate those re- 
lated to single-particle states (m*) with currently used 
functionals. This can be better understood by examin- 
ing the expressions for SNM and PNM EOS: 

+ C;|(^) 2/ V, (5a) 

E I T T\ 3 ^ f? 2\2/3 2/3 

+ [CZ(p ) + C?(po)] Po 

+ n+ci]\(^f%y\ (5b) 

If Cf(po) coefficients only contain one low power of the 
density (oc p\ ), the latter influences low-density parts 
of the EOS more than high-density ones. The effective 
mass term then determines the high-density part of the 
EOS. In SNM, this translates into the well-known rela- 
tion between m* and the incompressibility Koo 
In the case of PNM, the EOS above p sat is then mostly 
fixed by the term proportional to Cq + C{ in Eq. Il5b|) . 
and any attempt to use the density dependence to coun- 
teract its effects, results in a very strong constraint on 
the latter. This in turn degrades the behavior of the 
functional at and below saturation density and the fit to 
properties of finite nuclei. We recall at this point that the 
condition Am* > corresponds to C{ < 0, which drives 
the high-density PNM EOS down and explains why usual 
Skyrme functionals predict either a collapse of the PNM 
EOS if Am* > 0, or, like the SLy functionals fitted to 
PNM EOS, the wrong sign of the effective mass splitting 
in neutron rich matter. 

If Cf(po) coefficients contain an additional density de- 
pendence with a higher power, the previous discussion 
does not apply: using two density-dependent terms in 
the functional (cx p^ 3 ; p^ 3 ) ^3 allowed us to construct 
(/-, fo, /+) with a good fit to PNM EOS, a free choice 
of effective masses and satisfactory nuclear properties. 

The previous discussion already shows the type of 
problems and information arising from our attempt to 
improve on the fitting protocol of SLy functionals by us- 
ing more inputs from ab-initio calculations. Now, Fig. ^ 
shows SNM and PNM EOS as obtained from (/_, / , 
/+, SLy5) and as predicted by Variational Chain Sum- 
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mation (VCS) methods 51]. At this point, one can see 
that the four forces (/_, fo, f + , SLy5) reproduce both 
microscopic EOS with the same accuracy. However, it 
remains to be seen whether or not this translates into 
identical global spin-isospin properties and into similar 
nuclear structure properties. 




0.08 0.16 0.24 0.32 
P [fm 3 ] 



0.4 



1. Single-particle energies 

Effective masses are known to control the average den- 
sity of single-particle states. It is thus interesting to check 
to what extent such statement applies to neutron-rich nu- 
clei when varying m*. In this part of the study, we are 
mainly interested in evaluating the change in the single- 
particle energies generated by the functional for different 
splittings and not directly by a comparison with experi- 
mental results. 

Single-particle energies in 132 Sn and 208 Pb are plotted 
on Fig. |21 The general trend followed by neutron states 
with increasing Am* (from force /_ to / + ) corresponds 
to an increase of the density of neutron states: they tend 
to come closer to the Fermi energy ep; notable excep- 
tions being both neutron li levels in 208 Pb. The oppo- 
site behavior is observed in proton levels, which spread 
away from ep with increasing Am* (except for the proton 
l/in/2 level). However, these trends are rather marginal, 
which can be linked with the moderate bulk asymmetry 
of these nuclei (I = (N - Z)/A = 0.24 for 132 Sn and 
0.21 for 208 Pb). This moderate asymmetry means that 
the isovector term in the definition of the effective mass 
(Eq. Ijlbfl ) is weakly probed. 



Figure 1: SNM and PNM EOS as given by Skyrme functionals 
presently discussed (see text), compared with VCS results by 
Akmal et al. [HJ (crosses: PNM, plusses: SNM). 



D. Effects on properties of nuclei 

We now study the effects of the variation of the isovec- 
tor effective mass on selected properties of spherical nu- 
clei. We start with HF single-particle energies, then bind- 
ing energies, ending with a short sum-rule based analysis 
of isovector giant resonances. 

For computations of open-shell nuclei, we use, in the 
particle-particle channel, a zero-range interaction with a 
density dependent (mixed surface and volume) form fac- 
tor defined as: 



V pah (R,r) = V S(v) 



1 



pq(R) 

2p S at 



(G) 



where R = (i"i + f2)/2 and r = ri — T2- The local 
HFB equations are renormalized following the procedure 
developed in Ref. |H^ . 

The strength Vq is adjusted to the mean pairing gaps 
of six semi-magic nuclei (neutron gaps in 120 Sn, 198 Pb, 
212 Pb and proton gaps in 92 Mo, 144 Sm and 212 Rn). In 
this procedure we compute theoretical spectral gaps de- 
fined as A t h = Tr(fep)/Tr(/5), where h is the pairing field 
and p the pairing density |45j, and adjust each of them 
upon an experimental gap extracted through a five point 
difference formula from masses of neighboring nuclei, as 
suggested in Ref. [53| . 
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Figure 2: Single-particle energies [MeV] in 132 Sn and 208 Pb 
computed with indicated forces. 
Fermi level ef- 



Thick lines indicate the 



Let us therefore examine similar spectra for more 
neutron-rich nuclei, i.e. 78 Ni (/ = 0.28, experimentally 
observed [H) and 156 Sn (/ = 0.36). The 156 Sn nucleus is 
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used as an example of an extremely asymmetric system, 
even beyond the reach of planned radioactive beam facil- 
ities y^a]. We observe on the lower right panel of Fig.0 
that the effect of Am* on proton single-particle energies 
at Z = 50 is more pronounced in 156 Sn than it was in 
132 Sn. The modification of level densities appears quite 
clearly in 78 Ni also, while neutron levels around £f in 
156 Sn are shifted in a slightly more disordered way. 
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Figure 3: Same as Fig.^in 78 Ni and 156 Sn. 



High-£/low-n orbitals (n,£ being respectively the prin- 
cipal and orbital quantum numbers) are in fact more sen- 
sitive to variations of the spin-orbit field than to Am* be- 
cause of their spatial localization near the surface of the 
nucleus. The spin-orbit field is modified between func- 
tionals by the interplay between J 2 -term coefficients and 
effective mass parameters, since these both depend on 
the same non-local terms of the Skyrme force 56] . The 
spin-orbit force (pV ■ J terms in the EDF), which is sub- 
ject to a slight readjustment, does affect the spectra as 
well. We observed, overall, a marginal increase of the 
spin-orbit field strength when going from /_ to f + . This 
implies that while the global effect of modifying the level 
density is quite clearly observed when we alter the effec- 
tive mass parameters, details of the spectroscopy are at 
least as sensitive to the terms connected to the spin-orbit 
field. 



2. Pairing gaps 

As an example, neutron spectral gaps are plotted on 
Fig. 0] for Sn and Pb series, up to the drip line, against 
experimental gaps extracted through five-point mass for- 
mulas [4ll [53|. The slight change in the level den- 
sity translates into a modification of the pairing gaps: 
a higher neutron effective mass (/+) corresponds to a 
denser spectrum and higher gaps. The effect, which in- 
creases with asymmetry, remains however very small, be- 
cause of the limited alteration of single-particle levels 
seen on Figs. |21 and |3J 




10 20 30 40 50 60 70 80 90 100 
N-Z 

Figure 4: Neutron spectral gaps computed in Sn (bottom) 
and Pb (top) series with /_, fo, /+, as a function of asymme- 
try. Experimental A' 5 ' gaps extracted from masses |57j are 
plotted with error bars. 

In the end, the effect is negligible and would be over- 
whelmed by any other modification of the particle-hole 
part of the functional. For example the spin-orbit force, 
acting on the detailed level scheme, could alter the shape 
of gaps. The pairing functional itself is a subject of cur- 
rent debate regarding its density dependence, regulariza- 
tion scheme and finite-range corrections |58|, while the 
HFB formalism can itself be improved (particle number 
projection) as well as the mere choice of observable to 
be compared (definition of theoretical an experimental 
gaps) , although our choice has been proven to be possibly 
the most sound for extracting pure pairing effects |53| . 



3. Binding energies 

Let us now study the effect of the aforementioned vari- 
ation of level densities and pairing gaps on binding en- 
ergies. On Fig. we show the binding energy residuals 



Ui 



E exp for Sn and Pb isotopes and N = 50 and 
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N = 82 isotones. The evolution of E t h — E cxp along such 
series is usually plagued by an underbinding of open-shell 
nuclei with respect to closed-shell ones which translates 
into an arch shape of £7-residual curves. Although the 
variation of m* seems to impact the arches, again, the 
effect is negligible compared to the absolute value of de- 
viations from experiment, except in the N = 82 series 
where open-shell nuclei tend to be more underbound in 
the case of /+ . 
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Figure 5: Binding energy residuals computed with forces /_ 
fo and /+ for semi-magic series of nuclei, as indicated. 



4- Isovector giant resonances 

The isovector effective mass is usually defined from the 
energy- weighted sum rule mi (the TRK sum rule j4g) of 
the isovector giant dipole resonance (IVGDR): 



mi (£71 ; T= 1) = 



h 2 NZ 
2m A 



(1 + Kv) 



h 2 NZ m 
2m A m; 



(7) 



which exhibits its link with the strength distribution of 
isovector collective modes. We perform here a schematic 
study of dynamical pro pert ies of /_, /o, /+ by means of 
results derived in Ref. 59] . Thanks to RPA sum rules 
similar to Eq. Q , it is possible to fit an accurate param- 
eterization of the energy E\ — mi/m_i of isovector giant 
resonances in a given nucleus as a function of Skyrme pa- 
rameters. Results for GDR (L = 1) and isovector giant 
monopole (IVGMR, L = 0) modes in 208 Pb are shown 
in Table ITU compared to experimental energies (respec- 
tively from Refs. [6(| and j6J and corrected, as suggested 
in [53 , for the shift due to the spreading of the strength 
by damping effects - 2 MeV for GMR, 1 MeV for GDR). 



Table II: Ei energies of 208 Pb isovector giant resonances com- 
puted thanks to a sum-rule parameterization (see text), com- 
pared to experimental energy centroids. Experimental uncer- 
tainties are as indicated. We infer from figures in Ref. |5S^ 
the accuracy of theoretical energies computed with the fits in 
that reference, with respect to full RPA calculations, to be of 
the order of 1 MeV. 







£7i(L = 0,T = 


1) E 1 (L = 1,T=1) 


f- 


0.15 


24.55 


12.68 


fo 


0.43 


26.43 


13.60 


f+ 


0.60 


27.25 


14.01 


exp. 


centroid 


26.3 ± 1.1 


14.3 ± 0.1 



While /_ predicts both energies lower than experimen- 
tal ones, values for fo and /+ are compatible with exper- 
iment for the L = mode, and only / + approaches the 
experimental value for the L = 1 mode. This suggests 
that values of n v corresponding to a positive value of 
Am* (equal to or higher than 0.43 in our case) better 
describe isovector dynamics than lower values. 

As a summary, the effect of the splitting of neutron and 
proton effective masses with isospin asymmetry on single- 
particle energies, pairing gaps and binding energies, is 
noticeable and consistent, yet limited and thus hardly 
meaningful when compared to the overall (in)accuracy 
of the predictions made by current nuclear EDF. In fact, 
the main reason for not seeing a dramatic modification 
of EDF predictions when altering Am* is the limited 
amount of strongly asymmetric nuclear matter at high 
enough density in the ground state of nuclei with real- 
istic isospin as already suggested in [30j| . This makes 
the effect of the isovector effective mass rather marginal. 
Giant isovector resonances are certainly more fruitful to 
seek for an effect of a modification of Am*. Indeed, a 
sum-rule-based analysis of isovector collective modes al- 
lows a slightly more clear-cut conclusion, with a tendency 
to favor Am* > 0. The conclusion of the phenomenologi- 
cal study done in this section is that, while no observable 
listed here strongly ask for Am* > 0, there is no reason 
to omit this constraint in future functionals, since, as al- 
ready stated, ab-initio predictions for the sign of Am* 
are solid. There remains to check the intrinsic consis- 
tency of the functional in terms of other ab-initio inputs 
and stability criteria. 



III. FURTHER STUDY OF INFINITE MATTER 
A. Separation of the EOS into (S, T) channels 

1 . Introduction 

In this section, we discuss the contributions to the 
potential energy of SNM from the four two-body spin- 
isospin (S, T) channels. We compare our results with 
those predicted by BHF calculations using AV18 |6^| 
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two-body interaction and a three-body force constructed 
from meson exchange theory [&4l l65| . 

Using projectors on spin singlet and triplet states, re- 
spectively 



s=o 



5(1 -iU 



P 



S=l 



Pa 



(8) 



where P„ is the spin-exchange operator, and similar ex- 
pressions for isospin projectors Pt using the isospin ex- 
change operator P T , yields the potential energy in each 
(S, T) channel 
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kl 



PkkPlli 



(9) 



where the sum on k, I runs over all HF single-particle 
eigenstates whereas pkk designates the diagonal one- 
body density matrix. The notation \kl) denotes a non- 
normalized but antisymmetrized two-body state. In or- 
der to compare different many-body approaches (ab- 
initio or EDF), we use the "potential energy" which refers 
to the total binding energy from which is subtracted the 
kinetic energy of the non-interacting particle system. 

Note that due to the zero-range character of the 
Skyrme force, together with at most second-order deriva- 
tive terms, only L — 0, 1 partial waves occur explicitly 
whereas higher partial waves contribute to ab-initio EOS. 
We find, for SNM, 
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where (U,Xi) are usual coefficients of the Skyrme EDF, 
whereas (toi , xoi ) characterize the density-independent 
zero-range term and the two density-dependent ones, fol- 
lowing the parameterization of equation ijAlfl . 

The coefficients occurring in Eqs. I|10al) (|10d|l stem 
from the antisymmetrization condition (— ^ L + s + T = —1^ 
the relative angular momentum L being even for toi and 
t\ (k 2 ) terms and odd for t<z (k'-k) terms. The expression 
of the potential energy in channels (S, T) = (0, 0) and 
(1,1) is very simple since only the t 2 term contributes. 



2. Force vs. functional 

Previous statements, however, apply only to the case 
where the EDF is computed as the expectation value of 
an (antisymmetrized) effective force. In the more general 
case, it is still possible to define (S, T) channels starting 
from any Hartree-like functional. Indeed, the functional 
can always be expressed in terms of an effective non- 
antisymmctrizcd vertex and one can still plug a projector 
in the calculation of its matrix elements. In the pure 
functional case, there is however no more clear definition 
of partial waves, and spin-isospin channels emerge from 
the balance between coefficients of (iso)scalar/ (iso) vector 
couplings (see appendix[B]for a more detailed discussion) . 



As long as there are not enough inputs to constrain all 
degrees of freedom of a general functional, the "force" 
approach remains as an acceptable path, and hence shall 
be used in the following. 



3. Results 

Results are plotted against BHF predictions on Fig.[S] 
First, one can observe that results are rather scattered. 
Second, the main source of binding, from (S, T) — (0, 1) 
and (1, 0) channels, is not well described and the detailed 
saturation mechanism is not captured. It is clear that, 
even though all four forces reproduce perfectly PNM and 
SNM EOS, they do not have the same spin-isospin con- 
tent, and that the latter is in general rather poor. Thus, 
fitting the global EOS is an important element but it 
does not mean that spin-isospin properties of the func- 
tional are fixed once and for all. One needs to do more 
and fitting ab-initio predictions of Ep Q [ seems to be a 
good idea in the near future. However, one needs to make 
sure that the theoretical uncertainty of the data used is 
smaller than the expected accuracy of the fit to them. 
This calls for predictions from other ab-initio methods 
using the same two-body plus three-body Hamiltonian. 
Then, those ab-initio calculations should be repeated us- 
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ing different sets of two-body plus three-body Hamilto- 
nians in order to provide a theoretical error bar on those 
predictions. 
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Figure 6: Energy per particle in each (S, T) channel for SNM, 
as a function of density. Crosses refer to the BHF calcula- 
tions 621. 



The most obvious discrepancy appears in (0, 0) and 
(1,1) channels where Skyrme and BHF data have oppo- 
site signs above saturation density. The SLy5 parameter 
set shows a particular behavior in channel (1, 1) due to 
the choice of xi = — 1 to prevent ferromagnetic instabili- 
ties in PNM. Note that these two channels are taken care 
of, in the Skyrme functional, by the density-independent 
P-wave term only. The upper-right panel of Fig.0points 
out the tendency of Skyrme parameterizations to be at- 
tractive in polarized PNM, and hence to cause a collapse 
of its EOS at high density. At lower densities, BHF data 
show a distinctive behavior, being slightly attractive be- 
low p sat and repulsive above. This feature cannot be 
matched by the standard Skyrme functional which ex- 
hibits a monotonous behavior as a function of density in 
this channel, regardless of the value of (£2,^2)- 

It is also worth noticing that the failure in channel 
(1,1) becomes more and more prominent as one makes 
Am* closer to the ab-initio predictions (force / + ). The 
effective masses being governed by the momentum depen- 
dent terms of the force, it is not a surprise that the mod- 
ification of the former impacts channels (0, 0) and (1,1). 
What changes in the coefficients entering Eqs. (|10allUd|) 
stems only from the variation of m* and the associated 
rearrangement of parameters in the functional, most no- 
tably the C^f coefficients closely related to surface and 
surface-symmetry energies. The relatively tight require- 
ments on the latter imply that the four parameters of the 
non-local terms in the standard Skyrme energy functional 
would be dramatically overconstrained if we were to add 
the (S, T)-channel decomposition in the fitting data. 

In the end, the rather poor properties of the functional 
in channels (0, 0) and (1, 1), the degradation of the latter 
as the effective mass splitting is improved, the idea of us- 



ing ab-initio (S,T) contributions in the fit, call, at least, 
for a refinement of the odd-L term in the sense either of a 
density dependence or of a higher-order derivative term. 
The latter being prone to numerical instabilities and in- 
terpretation problems, a density-dependent k' • k term 
remains as one of the next potential enhancements to be 
brought to the Skyrme EDF (density-dependent deriva- 
tive terms have been considered already, but with a focus 
on even-L terms of the form ^(k 2 + k /2 )/?Q [66]). Phc- 
nomenological constraints on gradient terms are mainly 
related to the surface of nuclei, i.e. low-density regions. 
One can expect that, to first order, BHF data in channel 
{S,T) — (1, 1) can be matched with an extended func- 
tional while retaining a good agreement with other (ex- 
perimental) data. It is less clear in channel (0, 0) but 
further exploration of the extended parameter space may 
bring Skyrme and BHF data in better agreement. 



B. RPA linear response functions and the 
diagnosis of instabilities 

We attempt here to study general stability conditions 
of SNM with respect to finite-size density, spin, isospin 
and spin-isospin perturbations. Our basic ingredient is 
the RPA response function [67j derived analytically in 
Ref. j6^ for the central part of the Skyrme interac- 
tion. Recent work was done to incorporate the effect 
of the spin-orbit part which was found to be quite negli- 
gible |69|, and will be omitted in the present work. One 
starts by defining a one-body perturbing operator 



Q 



^ a 1 



(ID 



where a indexes particles in the system. The one-body 
spin-isospin operators Oa are defined as 
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1, 



e v a v =<r a T a , (12) 



where we use the denomination of (iso-)scalar (s) and 
(iso-)vector (v) channels in order to distinguish the un- 
coupled spin-isospin channels from the coupled two-body 
(S, T) channels discussed in the previous section. In 
Eq. (|El|) and the following, the first (second) subscripts 
denotes the spin (isospin). We then study the response 
to each type of perturbation separately through the re- 
sponse functions 



1 



1 



W - E n0 +ir) uj + E n0 - in 



(13) 



at the RPA level, where stands for a normalization 
volume and |n) is an excited state of the system, E n0 be- 
ing the corresponding energy. Since the central residual 
interaction does not couple the channels defined through 



10 



Eq. (|12fl in SNM, we can indeed consider each channel 
separately. 

The response function \^ can he seen as the propa- 
gator of the collective perturbation, i.e. the positions of 
its poles in the (q, lj) plane yield the dispersion relation 
of the mode. In this formalism, the onset of an unstable 
mode is marked by the occurrence of a pole in \ a ^ 
u> = 0, corresponding to zero excitation energy. Such a 
pole marks the transition between stable (x^ < 0) and 
unstable (x > 0) domains. Unstable modes of infi- 
nite wavelength (q = 0) arc those traditionally discussed 
in the context of Landau parameters. A pole at finite 
q characterizes a system which is unstable with respect 
to the appearance of a spatial oscillation of a given type 
(density, spin, isospin or spin-isospin) with a given wave- 
length A = 2n/q. 



The evaluation of response functions calls for the resid- 
ual interaction V p _h, defined as the second-order func- 
tional derivative of the energy with respect to the den- 
sity matrix. Its momentum-space matrix elements can 
be written, using total momentum conservation, as |6S| ]: 

V r p -h(qi,q 2 ,q) = (qi Q2 +q| Vp_ h |qi +q q 2 ), (14) 
= Wx^+Wiiq) (qi -q 2 ) 2 , (15) 

with 

W x {q) = i [ Wl s (q) + WF{q) cr 1 -a 2 + W^iq) n o f 2 
+ WT(g)o-i-<r 2 Ti°7H (16) 
and a similar expression for Wi. We find: 



Wf{q) 



Wr(q) 



Wr(q) 



wr(q) 



i(« + 6)(« + 3) i/3 A p 2 

n Po ~ 2 ^o q 



i=l 



2C- + 2^C-Vo /3 -2C As q 2 , 
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2Cf + 2 C?% /3 " 2Cf p q 2 , 
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= 2Cf + 2 Cx'pt - 2C As q 2 , 



(17a) 
(17b) 
(17c) 
(17d) 



Wi s (q) 



Wi v (q) 



With the above expression for the residual interaction, the response function reads as 



X {a) (w,q) =4n 



1 - W^Uq - 2W^ a 'k£ I q 



(a) 1.2 / t;2 



m'kl w (a) 
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Co, (18a) 

Cf, (18b) 

CT, (18c) 

Cf. (18d) 



(19) 



rn 



where g = q/2kp, v — ujm*/qk-p and 110,2,4 are general- 
ized Lindhard functions, see Ref. |68| . 

As already said, the limit q — > corresponds to per- 
turbations of infinite wavelength, keeping the system ho- 
mogeneous. There, the residual interaction is uniquely 
determined by Landau parameters Fi,F(,Gi,G[, with 
I = 0, 1, and well known stability conditions are obtained 
under the form 70]: 



1 



X, 



21 + 1 



>0, 



(20) 



where X\ represents any of the Landau parameters. We 
have used this criterion in the fit of our forces f x , ensur- 
ing that no spin or spin-isospin instability would occur 



below 2p sat . We observe that, from the point of view 
of Landau parameters, the most critical channel is the 
vector-isovector one, with associated instabilities at den- 
sities as low as 2p sat (see the upper-right panel of Fig.|5J). 
This behavior is linked to the attractive character of the 
functional in channel (S, T) = (1, 1) which gives rise to a 
collapse of spin-polarized PNM, and accordingly, a van- 
ishing spin-isospin symmetry energy. Therefore, better 
reproducing the decomposition into (S, T) channels of 
EOS obtained from ab-initio methods is not only a mat- 
ter of microscopic motivation, but also a necessity to 
avoid unwanted instabilities. 

Beyond infinite wavelength instabilities, we also aim at 
demonstrating that a more general treatment is needed 



11 



Table III: Values of the Cf " coefficient, in MeV fm 5 





f- 


SLy5 


fa 


/+ 


LNS SkP 


Am* 


-0.284 


-0.182 


0.001 


0.170 


0.227 0.418 


c a p 


5.4 


16.7 


21.4 


29.4 


33.75 35.0 



to fully describe and control unstable modes which arise 
in the Skyrme EDF framework. Thus, contributions to 
the residual interaction coming from functional terms of 
the form pAp are zero for q = 0, whereas such terms 
drive finite-size instabilities. 

Indeed, we have observed that existing (SkP) or new 
parameterizations built with a high value of k v in order 
to reproduce the microscopic splitting of effective masses, 
tend to spatially separate protons from neutrons in spher- 
ical mean-field calculations, where enough iterations lead 
to states with strongly oscillating densities and a diverg- 
ing energy. Following a preliminary phenomenological 
reasoning, we could relate this effect to the C^ p pi ■ Api 
term in the functional, as this term can energetically fa- 
vor strong oscillations of the isovector density pi which 
arise in the case of such a spatial n-p separation. More- 
over, Eqs. I|17all8d[) show that such a term can yield an 
attractive contribution to the residual interaction in the 
case of a short-wavelength (high q) perturbation. We 
found empirically that parameter sets for which this in- 



stability arises are characterized by a high value of C^ p ' , 
that is C^ p > 30. As seen from Table UTTl this parame- 
ter is strongly correlated with the effective mass splitting 
Am* in such a way that a positive splitting, as required 
by ab-initio predictions, leads to instabilities. 

Whereas we were obviously unable to provide a fully 
converged (and hence physically meaningful) and clearly 
unstable force to illustrate the previous statements, 
we found that certain forces available in the literature 
present the aforementioned behavior. For example, con- 
vergence problems have arisen (and have already been 
pointed out in another study ) for the SkP parameter 
set 0|- The nature of the instabilities discussed here 
is illustrated on the lower panels of Fig. where neu- 
tron and proton densities are plotted at various stages 
of execution of a self-consistent iterative procedure with 
SkP in 56 Ni. We see that strong, opposing oscillations 
of neutron and proton densities are formed, and steadily 
increase with iterations. Such a behavior happens after a 
seemingly converged situation for which the relative en- 
ergy variation is small but almost constant over a large 
number of iterations and the evolution of the energy is 
monotonous. 

The study of the linear response function in the scalar- 
isovector channel allows us to provide a more quantitative 
ground to the previous observation. By plotting criti- 
cal densities (lowest density of occurrence of a pole in 
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Figure 7: Neutron and proton densities in central regions of 
56 Ni (bottom) and 40 Ca (top) plotted for a fully converged 
computation using the SLy5 force (solid line; relative varia- 
tion of energy between iterations less than 10 -14 ) and along 
a series of iterations done with SkP (bottom) and LNS (top). 
Number of iterations as indicated in key. In both cases the col- 
lapse happens after a seemingly converged situation (~ 10~ 9 
relative energy variation, steady over a large number of it- 
erations indicating a nearly linear evolution of the energy), 
which can be mistaken for an energy minimum if too loose a 
convergence criterion is used. 



X^(uj = 0, q)) for a given q on Fig. EI we see that these 
critical densities can be very close to p sa t for q 2.5 to 
3 fm -1 . This is the case for SkP, which displays the low- 
est scalar-isovector critical density of all forces studied in 
this paper. Accordingly, it is the most prone to a lack of 
convergence in HF calculations. 




Figure 8: The lowest density of occurrence of a pole in 
X sv (t<J = 0,q) is plotted against the wave- number q of the 
scalar-isovector perturbation. The curves end at q — 2/cf 
since the ground state can not couple to excitations with 
cj = 0,q > 2k F . 



12 



The link between response functions and convergence 
problems can indeed be understood by classifying them 
by their magnitude: in case of a stable but very soft 
mode, lack of convergence arises from the existence of a 
continuum of quasi-degenerate mean-field states, among 
which no minimization or self-consistency algorithm shall 
be able to decisively find an energy minimum without a 
considerable amount of iterations. If the soft mode be- 
comes unstable, it causes a divergence of the energy and 
of other observable such as the densities. We see in the 
agreement between the RPA study of SNM and the obser- 
vation of unstable HF calculations of nuclei a qualitative 
validation of our Local Density Approximation (LDA)- 
based treatment of instabilities: soft or unstable modes 
occurring in INM happen for the same parameter sets in 
finite nuclei. 

The large number of iterations needed for the diver- 
gence to occur on Fig. is a consequence of the limiting 
case embodied by SkP, such that the existence of a def- 
inite instability is highly dependent on finite size effects 
(choice of the nucleus) and discretization details in the 
numerical procedure. If SkP is a limiting case, LNS also 
displays a low critical density in the scalar-isovector chan- 
nel (Fig. |8J. In this case, we observed proton-neutron 
separation in 40 Ca and for small mesh steps (0.1 fm) 
only (see Fig. 01 , while it is more frequent with SkP. Our 
force /+, with a critical density just slightly higher than 
LNS, successfully passed the test of computing a series 
of 134 spherical nuclei. This again demonstrates that 
testing finite-size instabilities through response functions 
constitutes an accurate tool, the critical density (and its 
proximity to p sa t) being a good measure of the gravity 
of the problems one might encounter in finite nuclei. Al- 
though the actual occurrence of instabilities is subject to 
details of the numerical treatment, it is now clear that 
their origin can be traced back to the choice of parame- 
ters in the functional itself. 

Nevertheless, even if a functional does not display clear 
instabilities but only spurious soft collective modes, con- 
vergence difficulties shall arise in mean-field calculations 
while such a mode will translate into a non-physical 
low-lying spectrum in a beyond-mean-field framework. 
This can then yield excessive correlation energies if one 
systematically includes correlations in the ground state 
e.g. in (Q)RPA or GCM methods. One should thus make 
sure that no spurious (even remotely) soft mode occurs 
at saturation density in order to prevent such problems. 

Having demonstrated the importance of finite-size in- 
stabilities, let us go back to discussing our original set of 
forces and perform a generalization to other spin-isospin 
channels. 

Critical densities are plotted on Fig. [5] for the four 
channels defined in Eq. Q12[l . The upper- left panel shows 
that, while no unstable mode occurs at q = thanks 
to fitting PNM EOS to relatively high density, scalar- 
isovector instabilities may happen little above p sa t for 
q w 2.5 to 3 fin - . In addition, there is a clear trend 
for lowering the critical density when Am* is increased, 



in agreement with the preliminary phenomenological rea- 
soning on C^ p . 
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Figure 9: Same as Fig. |H| for all spin-isospin channels. The 
lower-left panel shows the region of spinodal instabilities be- 
low psat. The domain of q covered in this case determines the 
size of structures formed, while the region between 0.1 and 
0.16 fm -1 appears as metastable. 

Spin channels have been taken care of during the fit 
thanks to Landau parameters, which describe the resid- 
ual interaction at q = 0. The result can be seen on the 
right panels of Fig. [§1 where the critical densities of in- 
stability are plotted for spin-flip modes (isoscalar and 
isovector). As previously stated, the most dangerous 
q = instability is found in the vector-isovector chan- 
nel. By looking at the upper-right panel of Fig. [5] one 
can see that the critical density is even lower (1.5p sa t) at 
q = 3 fm - than at q — 0, a domain not covered by the 
criterion of Eq. (J2DJ. 

An even more prominent finite-size effect can be ob- 
served in the isoscalar spin-flip channel (lower-right panel 
of Fig. where, while no instability occurs at q = as 
in the case of most Skyrme forces, finite-size instabilities 
occur at low density, even below p sat ! These instabil- 
ities are linked to the C$ s So • Aso term which makes 
the vector-isoscalar V^-h attractive at large q whereas 
it is repulsive at q = 0. Values of C^ s , indeed, are as 
high as 45.85 and 47.32 for SLy5 and /_, respectively. 
As a consequence, one can expect divergences in calcu- 
lations of odd or rotating nuclei with the latter forces 
if the aforementioned terms are included. In this case, 
though, increasing Am* pushes the critical density far- 
ther from p sa t: fo and /+ functionals are thus the only 
ones to be free from instabilities at p sa t, without being 
totally satisfactory either. 

The previous discussion is valid if the full time-odd 
functional is taken into account. This must be stressed 
since So • Aso terms, which drive the most critical, 
finite-size instabilities, have never been included in self- 
consistent mean field calculations to date. However, 
RPA calculations are commonly performed by computing 
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the residual interaction matrices directly from the an- 
tisymmetrized force (plus rearrangement terms), which 
amounts to implicitly including the contribution to V^-h 
from all terms in the functional 

The latter findings finalize the picture of a competi- 
tion between spin and isospin instabilities. All in all, 
the strong interplay between the various quantities linked 
to the four parameters of the non-local terms in the 
Skyrme force does not seem to allow for a fully satis- 
factory compromise between stability criteria and ab- 
initio constraints on Am*. Again, we see that the non- 
local part of the Skyrme force is too simplistic to control 
all relevant properties. An extension with density- and 
momentum-dependent terms, allowing the fine-tuning of 
the functional at various densities, combined with the 
formal checks advocated in this paper, could prove to sig- 
nificantly improve the predictive power of Skyrme EDF. 



IV. CONCLUSION 

We have built a series of Skyrme energy density func- 
tionals to study the effect of a variation of the split- 
ting of neutron and proton effective masses with isospin 
asymmetry on properties of this mean-field-based model. 
Thanks to the use of a second density-dependent term 
in the underlying effective force, we could cover a wide 
range of effective mass splittings (Am*) with a satisfac- 
tory fit to nuclear properties. Indeed, nuclear observable 
predicted by our functionals /_, /o and /+ show a re- 
markable similarity, pointing out that spectra, pairing 
gaps and masses of bound nuclei are weakly sensitive to 
Am* , mostly due to their relatively low isospin asymme- 
try. Although observable were affected in a noticeable 
and consistent way, no clear improvement was seen when 
altering Am* either way. 

Beyond this phenomenological study, we have com- 
pared the splitting of the equation of state of symmet- 
ric infinite matter into spin-isospin channels provided by 
our functionals and by ab-initio Brueckner-Hartree-Fock 
calculations. Such a comparison showed an obvious dis- 
crepancy in (S,T) = (0,0) and (1,1) channels, where 
energies predicted by Skyrme functionals and by BHF 
calculations have opposite signs. The inconsistency in 
channel (5, T) = (1,1), where the Skyrme functional is 
attractive, translates into a collapse of polarized neutron 
matter EOS, related to the onset of spin-isospin instabili- 
ties at quite low density (2p sat ). In this channel, ab-initio 
predictions cannot be matched (in the Skyrme "force" 
approach) without an extension of the P-wave term. We 
also identified finite-size isospin instabilities caused by 
strong isovector gradient terms, which prevent the con- 
vergence of mean-field calculations. We were able to pro- 
vide a firm and quantitative basis to these observations 
through an analysis of finite-size instabilities by use of 
RPA linear response functions in SNM. The latter showed 
that finite-size effects in the analysis of instabilities tend 
to always dominate. 



The present study leads us to propose the systematic 
inclusion of consistency checks with ab-initio predictions 
of spin-isospin properties in the construction of our future 
functionals, as well as a systematic diagnosis of finite-size 
instabilities. 
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Appendix A: SKYRME ENERGY FUNCTIONAL 

We take the particle-hole part of the functional as given 
by the expectation value of a Skyrme effective force in- 
cluding two density-dependent terms: 

2 

V(R, r) = £ oi (l + xoiP*) S(r) [p (R)f 3 

1=0 

+ i*i(l + a:iP a ) [S(r) k 2 + k' 2 S(r)} 



+ t 2 (l + x 2 Pa) k' • *(r) k 

+ iW [<ri + <t 2 ] k' x 6(r) k, (Al) 

with the usual notations 

R = (n+r 2 )/2, (A2a) 

r = n-ra, (A2b) 

k = h (Vi_ V2) ' (A2c) 

k' = C.C. of k acting on the left, (A2d) 

er = o-i +er 2 , (A2e) 

Pa = \ (1 + ffi-o-a). (A2f) 

The total binding energy of a nuclear system can be 
written as a functional of a local energy density 

E = J d 3 r H(t), (A3) 

h 2 

H = - — To + Wskyrme + Wcoul., (A4) 



U -7/cvcn i -i/cvcn , i_/odd i oyodd / a r\ 

AtSkyrme — tt + ri 1 + rt + rt 1 , (AO) 

where the superscripts in the last equation indicate the 
behavior with respect to time reversal of densities occur- 
ring in each term, while subscripts indicate the rank of 
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the densities in isospin space. The corresponding expres- sions are 
I 



ni ven = OS pi + C* p Po A Pa + Cl Po t + C J J 2 + C VJ p V ■ J , (A6a) 

HT cn = G{ p? + Cf p pi ■ Ap! + Cl p x ■ ri + C( 1^ + Cj 3 pi-V-Jx, (A6b) 

H° dd = C S s 2 + C As s • As + Cf s • To + C Vs (V • s ) 2 + C 3 Q j 2 + s • (V x j ), (A6c) 

H°x dd = C{ + C As s x • Asi + Cf si • T x + C^ s (V • s^) 2 + C{ j x a + Cf J " s r (Vx ji). (A6d) 



Bold letters denote vector densities and arrows denote 
isovector densities. Neutron and proton densities (q = 
+1 for neutrons, q = — 1 for protons) are thus given by 



Pq 0) 
r q (r) 
l q (r) 



(po + qpx,a), 
(r + g n, 3 ) , 
■ (Jo + 9 Jl,3) , 



(A7) 



with similar expressions for time-odd densities. The spin- 
current vector J t is built from the antisymmetric part of 
tensor J t . 

Let us give the expressions, in terms of Skyrme force 
parameters, of the coupling constants entering the HFB 
calculations which are altered by the addition of a second 
density-dependent term 



r<i> 



C{ 



e c o i Po 3 - e § toi 

E c p / P f 

i=0 
2 1 



/3 



(A8a) 



(A8b) 



i=0 



as well as the constants related, respectively, to isoscalar 
and isovector effective masses, 



1 

16 
1 

16 



[3ii + t 2 (5 + 4a:2) ], 



(A9a) 



[-tx(l-2xx)+t 2 (l + 2x 2 )], (A9b) 



and constants multiplying gradient terms discussed in 



sec. Illl Bl 
















-I, 

64 1 




+ t 2 {5- 


bix 2 ) ] 




(AlOa) 


C Ap = 


1| 

64 1 




+ 2xx) 


+ h(i- 


^2x 2 ) 


],(A10b) 




1| 
64 


3tx(l 


-2xx) 


+ t 2 (l- 


\-2x 2 ) 


] , (AlOc) 


C* As = 


1| 
64 


; stx + 


t 2 }- 






(AlOd) 



The expressions of all other coupling constants are given 
in Ref. [lj. Some of the above constants are linked 
through the gauge invariancc of the functional, related to 
the Galilean invariance of the underlying effective force: 



c{ = -ci, c t J = -c? T , c7 3 = c? J . (ah) 

The densities used above can be written as functionals 
of the density matrix expressed in coordinate space, i.e. 

p(xaq,x'a'q') = ^{k\x' a' q')(xaq\k) p kk , (A12) 



as 



j(r) = J d 3 xd 3 x' E S (* - x ) S (*' - x ) <V« <W 

aa'qq' 

x p(xaq,x'a'q'), (Al3a) 
Apo(r) = y"d 3 xd 3 x' E S(r-x)S(x' -x) <V, <W 

x (V' 2 +2V'- V + V 2 )p(xaq,x'a'q'), (A13b) 
r (r) = y"d 3 xd 3 x' E S{t-x)6(x' -x) 5 q - q <W 

<t cr' qq' 

x V' • V p(xaq,x'a'q'), (Al3c) 
Jo(r) = y"d 3 xd 3 x' E <5(r - x) <5(x' - x) fy, 

x rr(V- V')(8o- CT v p{xaq,x'a'q'), (A13d) 
jo(r) = y"d 3 xd 3 x' E <5(r - x) 5(x' - x) <V 8 <W 

(A13e) 



aa qq 



x — (V — V ) p(xaq,x'a'q' 
2i 



where cr, cr' are indices referring to spin, q and c/' refer 
to isospin, V is the gradient operator acting on the co- 
ordinate x, V' being the same acting on x' . Isovector 
and other time-odd densities can be expressed by replac- 
ing, respectively, S q > q by f q > q and S a > a by <r a > a where 
appropriate. 
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Appendix B: SEPARATION OF THE ENERGY 
INTO SPIN-ISOSPIN CHANNELS 



When the EDF is defined as the expectation value of 
an effective Hamiltonian, separating it into spin-isospin 
channels is straightforward, as in Eq. However, one 
can extend this definition to the case of any Hartree-like 
functional: let us start by recalling that in the case of 
the Skyrme force, the direct and exchange terms have 
the same analytical structure; one thus usually uses the 
expressions 



E pot = r ^ ( kl |Vskyrmc| M ) p k k Pll, 



(Bl) 



kl 



\kl) = \kl)-\lk) = (1 - P r P a Pr) \kl) , (B2) 

where the last expression uses the position, spin and 
isospin exchange operators to define an antisymmetrized 
and non-normalized two-body state. One then writes 
down the antisymmetrized form of Eq. IjAljl and the EDF 



by using the definition of densities entering Eqs. (|A6ajl - 

Leaving the antisymmetrized Hamiltonian framework, 
it is always possible to define the potential part of the 
functional as the direct term of the expectation value of 
a certain operator, as in 



£ P ot = ^2(kl\V EDF \kl 



Pkk Pit, 



(B3) 



kl 



recalling that Pedf = V^kyrme(l — PrPa Pt) hi the Hamil- 
tonian case. One then defines the energy per channel as 



E 



ST 
EDF 



Vedf PsPt 



kl ) pkk Pll, (B4) 



which, with the definitions (|A6a|) - (|A6d|) for coupling con- 
stants, yields (retaining only terms acting in infinite mat- 
ter) 



Eg = J d 3 r TL ST (r) (B5) 

H ST = [C P + (45 - 3)Q + (AT - 3)Cf + (45 - 3)(4T - 3)Cf] 

x -L [(25 + 1)(2T + l)p 2 + (25 - 1)(2T + ljs 2 , + (25 + 1)(2T - l)^ 2 + (25 - 1)(2T - l)^ 2 ] 
16 

+ [C T + (45 - 3)C sT + (AT - 3)C{ + (AS - 3) (AT - 3)Cf] 

x -^[(25+l)(2T+l)p To + (25-l)(2T+l)s -T + (25 + 1)(2T - + (25 - 1)(2T - l)sr • f x ] . 

16 

(B6) 



Appendix C: PARTICLE-HOLE RESIDUAL 
INTERACTION 



With the definition of densities given in Eqs. lA13b|) - 
(|A13e|) the particle-hole residual interaction is obtained 
through 

<t[ q[ , x 2 <t' 2 q' 2 \ V p - h | x x a x qy , x 2 cr 2 <?2 ) 
S 2 £ 



5p(yi 1 a 1 q 1 ,x.' 1 a' 1 q' 1 ) <5p(x 2 cr 2 q 2 , x 2 cr 2 <j 2 ) ' 

which, for the central, spin-scalar, isoscalar part of the 
functional, i. e. 



C P M Pt + C* p p A po + C T (po T - jg), 



(C2) 



(the generalization to the full central part, omitted here 
for the sake of brevity, being straightforward), reads 



v^ h = ics* + [{ + 2 U + 1 )p° 



-G 



-2(Vi-v 1 )-(v 2 -v 2 ; 

C Ap ((V' 1 2 + 2V' 1 • Vi + V 2 ) 



-(V 2 2 + 2V 2 • v 2 + v 2 ; 



i/3 



(C3) 



When computing the momentum-space matrix element, 
Eq. I|14f> . one makes the substitutions Vi = ?'(qi + q), 
V 2 = iq 2 , Vi = -iqi an d V 2 =_£(q 2 +q) (with H = 1), 
which yields expressions (|17a|) - (|18d|) . 
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